function [beta,Var,resid] = nwest_iv(y,X,Z,lags)   
[nObs, nVar] = size(X); 

%% Point estimate: drop missings
notmissing = all([~isnan(y) ~isnan(X) ~isnan(Z)],2); 
Xnm = X(notmissing,:); 
Znm = Z(notmissing,:);
ynm = y(notmissing,:);

beta      = (Znm'*Xnm)\(Znm'*ynm);
yhatnm    = Xnm * beta;
residnm   = ynm - yhatnm;

%% Standard errors: restore correct time structure, assume unobserved moments are zero
resid = zeros(nObs,1);
resid(notmissing) = residnm;

Zm = Z;
Zm(isnan(Z)) = 0;

E = [];
for i=1:nVar
    E = [E; resid'];
end
       
ezzpm = E.*Zm';
V = zeros(nVar,nVar);   
V = V + ezzpm*ezzpm';

for h=1:lags+1
    vh = zeros(nVar,nVar); 
    Gh = ezzpm(:,(h+1):nObs) * ezzpm(:,1:nObs-h)';
    vh = vh + Gh + Gh';
    V = V + (1 - h/(lags+1)) * vh; 
end
    
Var = (Znm'*Xnm)\ V / (Znm'*Xnm);